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Abstract 

In this paper two-dimensional electromagnetic scattering problems with a time-periodic 
incident held are considered. The scatterer is a perfect conductor, and an artificial boundary 
condition is used. The large time behavior of solutions, depending on (divergence-free) 
initial conditions, is characterized. It turns out that in addition to the expected time- 
periodic solution the limiting solution may also contain a spurious stationary held. The 
source of the stationary held is explained and equations describing it are obtained. Several 
avoidance strategies are discussed, and numerical comparisons of these techiques are given. 


‘This work was supported by the United States Air Force Office of Scientific Research under grant 
F-49620-94-f-03ff. Research was also supported by the National Aeronautics and Space Administration 
under NASA Contract No. NAS1-19480 while the authors were in residence at the Institute for Computer 
Applications in Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23681- 
0001 . 


1 



1. INTRODUCTION 


For solving scattering problems with a time-periodic incident field, the periodic solution is 
often obtained by choosing initial conditions and time marching Maxwell’s curl equations 
to a periodic state. The initial conditions must satisfy the divergence equations but are 
otherwise arbitrary. 

When this time marching method is used to solve scattering problems for perfect conductors 
it is found that the computed periodic solutions are “contaminated” by spurious stationary 
components. For example, in two dimensional transverse magnetic calculations it is observed 
that the magnetic field is affected by a spurious component; in three dimensions both the 
magnetic and electric fields are affected. The spurious fields can be relatively large and must 
usually be removed empirically by postprocessing, for example by peak-to-peak averaging 
[2], [3], or by the postprocessing method mentioned in section 4 below. These methods are 
usually effective and permit quantities such as radar cross sections to be obtained from the 
resulting numerical data. However, postprocessing introduces additional uncertainties and 
it is worthwhile to try to obtain accurate numerical solutions directly. 

In this paper we will explain the origin of the spurious fields and obtain the differential 
equations and boundary conditions which characterize them. Based on these results we will 
propose some techniques which enable accurate results to be found directly and efficiently 
without postprocessing. We will provide numerical comparisons of the techniques. 

The next section illustrates the occurrence of a spurious field in the relatively simple setting 
of a two dimensional transverse magnetic computation. Then we explain the source of this 
field and derive the equations which define it. Following that we discuss some avoidance 
strategies and illustrate their properties. The remainder of the paper extends the results to 
other problems and briefly mentions some mathematical issues which are analysed in detail 
in [4], 


2. THE MAGNETIC OFFSET 


In a two dimensional transverse magnetic problem it is the magnetic field which is affected 
by the spurious stationary component. To illustrate the spurious field we consider the 
following transverse magnetic scattering problem: in 0, 
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Figure 1: The time history of H a . (solid hue) and H y (dashed line). 


and initial conditions 

E(x, 0) = E 0 (x), H(x, 0) = H 0 (x). (6) 

In these equations, E and H denote the scattered held variables, E; is the incident held, 
r s denotes the boundary of the (perfect conductor) scatterer and r o is the artificial outer 
boundary, denotes the domain between the inner and outer boundaries, n is the outer 
normal to and x = (x,y). We assume that E; can be written in the form 

E i '(x,/) = e ^E i '(x), 

and we look for the time-periodic solutions of ( 1 )— ( 6 ) in a similar form. For simplicity, we 
are assuming a first order radiation condition on the outer boundary . The results are not 
materially affected if higher order radiation conditions are used. 

To illustrate the magnetic offset we take a square scatterer with sides aligned with the 
coordinate axes and electric size 1 based on half the side length of the square. The scatterer 
is illuminated from the (1,0) direction by a held of the form 1000e i ^ t_k ' x ^. The outer 
boundary is ^ wavelengths (4.5 body lengths) from the center of the scatterer. This 
problem was marched to a periodic state, starting from zero initial conditions, using the 
standard uniform mesh “FDTD” scheme [1] with 207T points per wavelength and with a 
timestep of l/\/2 the maximum size allowed by the stability condition. 

In Figure 1 we show a time history of H a . and H y sampled at the points ( — y, |^~) and 
( — 20 ^’ It )’ respectively (measured from the center of the body in wavelengths). While the 
records appear to approach periodic functions, they are not symmetric about the horizontal 
axis and consequently they cannot be of the desired form e iL ^ t H(x). A similar situation 
holds at the other mesh points, the amount of the offset being a function of position. This 
function is spurious and is not part of the solution to the periodically forced scattering 
problem. We wish to characterize it with a view to eliminating it altogether. 

In contrast to Figure 1, Figure 2 shows the time history of the electric held at the point 
( — y, y). There is no obvious problem with this solution and in fact it is close to the exact 
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solution of the periodically forced difference equations. A general proof of this is given in [4]. 


3. ORIGIN OF THE MAGNETIC OFFSET 


Let C^E and denote the periodic solutions obtained by solving the periodically forced 

differential equations without the initial conditions. We shall assume that the solutions E 
and H of the initial value problem are such that for large times, 

E(x,/) ~ e^E(x), H(x,/)~e^H(x) + H*(x). (7) 


These assumptions reflect wliat is seen computationally. 

Since C^E and satisfy the same equations as E and H (excluding the initial condi- 

tions) we obtain by subtraction the following equations for H*: 


div H* = 0 

curlH* = 0 

DU* 


H* x n|r D = 0. 


( 8 ) 

(9) 

( 10 ) 

(ID 


To obtain a boundary condition on the scatterer we integrate Faraday’s law with respect 
to time, use the initial conditions and evaluate the normal components on T s : 


H(x) • n = Ho(x) • n / curlE(x, r ) • n dr 

H Jo 


on T s . 


Since curlE • n on T s is the tangential derivative of E along the boundary, it is completely 
determined by the boundary condition. So 


H(x) • n = Ho(x) • n -| / curl E;(x, r ) • n dr 

/' Ji) 
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f 1 ~ \ e tujt 

= Ho(x) • n — curlEj-(x) • n + curlEAx) • n on T s . 

V tuy. J tuy. 

In the above expression it is only the second term which has the required periodic behavior 
and it follows that the first term, which is stationary, must correspond to the normal 
boundary condition for H*. Thus we are led to 

H* • n| r = (H 0 • n - — curlE, • n)| r • (12) 

s VjJfi 

To interpret this boundary condition, assume that the incident magnetic held is also of the 
form e* lx ' t Hj'(x) such that Maxwell’s equations are satisfied. In that case Faraday’s law gives 

e“ 4 H,(x) = — e’^curlE,- 

lUfl 

and evaluating this at / = 0 and using it to replace the second term in the boundary con- 
dition it follows that H* • n is the normal component of the total magnetic held at / = 0. 
Since B • n must be continuous across T s it follows that there is a magnetic held inside 
the scatterer if the right side of (12) is nonzero. In addition this held is stationary since 
E is zero inside a perfect conductor, and we may attribute it to a steady current howing 
unhindered in the conductor. It is this current which produces the magnetic offset. 


4. AVOIDING THE MAGNETIC OFFSET 

In this section we will discuss four different possibilities for obtaining the correct magnetic 
held e* lx ' t H(x). The hrst approach is a simple postprocessing technique which does not 
use the results of the previous sections except for the assumptions (7). The other three 
approaches do rely on the equations obtained in the previous sections. 

Method 1: Postprocessing. The postprocessing approach consists simply of htting a 
function of the form A + Be tujt to the computed H at each mesh point, and then discarding 
the stationary part A. A decision has to be made about when to sample H. This should be 
at reasonably spaced times after the solution is believed to have become adequately periodic 
(3-6 periods is often sufficient). 

To illustrate the stationary part we took a square scatterer as before and illuminated it 
from the (1,0) direction with a held E 4 -(x,/) = lOOOe*^ - ^). The problem was solved (for 
the imaginary part of the solution) by the method mentioned above in section 2, starting 
from zero initial conditions. We calculated (the imaginary part of ) the stationary part by 
sampling the values of H after m, m — j and m — j- timesteps, where T is the period of 
the incident held in timesteps, for m = 500 and m = 1000 (about 4 and 8 periods). The 
values of the stationary parts of the x- and r/- components of H on lines y = and y = 
respectively, are shown in hgures 3 and 4. Also shown are the theoretically predicted values 
of the stationary held H* (the approximate solution of (8)-(12) ). We see that the values of 
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Figure 3: The stationary part of H a .. 


Figure 4: The stationary part of H y . 


Postprocessing method (zero initial data) 

500 timesteps — solid line, 1000 timesteps — dashed line, H* — dotted line. 


H* are quite near to the stationary parts of H, especially for longer time. The stationary 
part is relatively large compared to the amplitude of H, which is near one for most of fi . 

Method 2: Solving the div-curl system. The second approach to elimination of the 
magnetic offset is based on the observation from (12) that if the initial condition for H 

satisfies Ho • n = curl E; • n on T s , then the div-curl system for H* will have only the 

/„•// 

zero solution. One way to choose a suitable Hq is as a solution of the equations 


curl H 0 = 0 (13) 

div H 0 = 0 (14) 

Ho • n = curl E; • n on T s (15) 

/„•// 

Ho - n 0 on 1,,. (16) 


To illustrate this we solved Maxwell equations with the initial conditions Eo = 0 and Ho 
being the solution of ( 14 )— ( 16 ). Then we calculated the stationary part as above. The 
results (the values of the stationary parts of the x- and ^-components of H on lines y = 
and y = y, as before) are shown in Figures 5 and 6. 

One can see that the stationary part is several times smaller than in the postprocessing 
method; in fact it is probably caused by the discretization error in time. However, even 
though there is a natural way to formulate the div-curl system in the FDTD framework, 
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Figure 5: The stationary part of H a .. 


Figure 6: The stationary part of H y . 


H 0 is a solution of the div-curl problem ( 14 )— ( 16 ). 

500 timesteps — solid line, 1000 timesteps — dashed line 


solving it can be quite expensive. Also, up to round-off error the solution is similar to that 
of the postprocessing method after discarding the stationary part. So it seems that this 
method is not worth while. 


Method 3: Extension of E;. Another way to satisfy the condition Ho-n = curlE 8 n 

iwf.1 


on r s is to use the initial conditions 


E 0 (x) = — ^(x)Ej-(x), H 0 (x) = T^curl(^(x)Ei(x)), (17) 

where cj) is a smooth real- valued function such that cj) = 1 in a neighborhood of T s , and cj) = 0 
in a neighborhood of r o . A possible advantage of this scheme is that since the initial values 
and the boundary data match at / = 0, the solution should be smoother in time for small 
times. We would then expect that discretizations should have a better rate of convergence 
as the mesh size approaches zero. 

We used for cj) a product of continuously differentiable cubic splines in x and in y. The 
resulting “stationary held” is shown on Figures 7 and 8. We see that it is much smaller than 
in the previous case, especially for longer time (40-80 times smaller after 1000 timesteps). 

Method 4: Modified incident field. This method consists of giving zero initial data 
and making an adjustment to the incident field. It is possible to change this field for a 
small interval of time so that the spurious field will be zero. One possibility is to define the 
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Figure 7: The stationary part of H a 


Figure 8: The stationary part of H^ 


The smooth initial data (17). 

500 timesteps — solid line, 1000 timesteps — dashed line 


incident held to be 


E;( x, / ) = 


_ J 2vr 


— e^Edx), 0 < / < — , 


e lujt E ,-(x), / > — . 

UJ 


Then at time / 0 = — , on T s we have 


H (x,/ 0 ) • n = - [ t0 —e iujt dt curl E,-(x) • n = — e^curl E,-(x) • n, (19) 
j-i Jo t 0 

so that H(x,/o) satisfies the correct condition at time to = yf- Again the initial and the 
boundary conditions match, so the solution is smoother. Since the incident field reaches the 
correct value after the first period, the additional work required is roughly equal to solving 
the equations during one period. 

The results are shown in Figures 9 and 10. 

For an incident field in the form e *M-k-x) ^ m jg^ b e preferable to adjust the field as 
follows: 

10, u>t < k • x, 

E,-(x,/)=J Ut ~ k ' x e «M- k - x) , 0 < ojt - k -x < 2tt, (20) 


2vr 


2ir < ujt. — k • x. 
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Figure 9: The stationary part of H a .. Figure 10: The stationary part of H y . 

Continuous adjustment of the incident held, given by (18). 

500 timesteps — solid line, 1000 timesteps — dashed line 


Again at time to large enough so that the incident held has reached the correct value on all 
r s , we get that H(x,/o) satisfies the condition necessary to make the spurious held vanish. 
This adjustment makes the imaginary part of E; continuously differentiable, so we should 
get even better convergence when solving for the imaginary part only. The results are shown 
in Figures 11 and 12. (We used m = 600 instead of m = 500, because after 4 periods the 
adjusted incident held has barely reached its correct value near the outer boundary.) Note 
that the stationary parts are the smallest compared with the other methods. 

To compare the smoothness of the solutions obtained by the different methods we present 
here the time history of H a . at the point ( — y, ). This point, as any point near the .r-axis, 
is special in the sense that the amplitude of H a . is very small, so effects of nonsmoothness 
stand out more clearly. Figure 13 shows the solutions with zero initial data (method 1) 
and in the case of “correct” smooth initial data (17) (method 3). The graph of the solution 
obtained by first solving the div-curl system (method 2) looks the same as of the solution 
obtained by method 1, only shifted by a constant (the size of the shift is determined by 
the value of the solution of the div-curl problem at given point), so we chose not to present 
it here. The solutions exhibit very large oscillations in the beginning. They correspond 
to a large error in approximating the discontinuous wavefront. As the front moves away, 
the solution becames smoother. Then the oscillations reappear approximately at the time 
when the (reflected) wavefront arrives back from the corners of the artificial boundary (the 
reflections from the smooth part of the boundary are probably too small to be seen). Figure 
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Figure 11: The stationary part of H a .. 


Figure 12: The stationary part of H y . 


Smooth adjustment of the incident held, given by (20). 
600 timesteps — solid line, 1000 timesteps — dashed line 


14 shows the solutions with zero initial data and altered boundary conditions (method 4, 
boundary conditions (18) and (20) ). We see that the solution looks much smoother than 
the solutions obtained with the other methods, so one might expect to have smaller dis- 
cretization error than in the other cases. 


9 




Figure 13: The time history of H a . at the point ( — in the case of zero initial data 

(solid line) and the smooth initial data (17) (dashed line). 



Figure 14: The time history of H a . at the point ( — y,^y) in the case of continuous ad- 
justment of the boundary conditions, given by (18) (solid line), and in the case of smooth 
adjustment (20) (dashed line), both with zero initial data. 


5. OTHER PROBLEMS 


Consider the transverse electric problem: in f 1, 


8E 

e ~di 

du 

div E 


curlH 

— curl E 

0 


with boundary conditions 


E x n = — E; x n on T s 

E x n + cf.iH = 0 on r o 
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and initial conditions 


E(x, 0) = E 0 (x), H(x, 0) = H 0 (x). 

In this case there is an electric offset. The magnetic offset is zero because H • n = 0 on T s 
holds automatically. 

Assuming that for large times 

E(x, t) ~ e^E(x) + E*(x), H(x, t) ~ e^H(x), 
we obtain as before by subtraction the following equations for E*: 

<9E* 

~df ~ ° 

curl E* = 0 

div E* = 0 

E* x n = 0 on T 0 and on T s 

A system such as this is determined up to a one parameter family of solutions. In this case 
the parameter is the total surface charge on the body. This quantity should be zero to avoid 
a spurious addition to the electric held. It is computed as follows: 

[ E • n = [ Eo-nH — [ [ curl H • n = [ Eo • n. 

.'i .'i < Jo J\ J\ 

It fohows that the total charge is constant in time and it is this charge, determined by 
the initial condition, which is responsible for the electric offset held. The simplest way to 
eliminate the electric offset is to put Eo = 0. The system above, supplemented with this 
condition of zero surface charge then has only the zero solution. 

For the three dimensional case the calculations above remain valid so both the electric and 
magnetic helds may have spurious components. Our equations for E* and H* continue to 
hold and so should the remedies although we have not checked that explicitly. 

It has not been possible to completely rule out the possibility of additional (possibly time 
dependent) spurious helds in three dimensions although we do not expect that such helds 
exist. 


6. MATHEMATICAL COMMENTS 

We have made a number of assumptions in order to obtain the equations for the offset helds. 
We wish to mention the major assumptions here. 

First there is the question of whether the equations (l)-(6) admit a periodic solution. This 
is equivalent to the question of existence for Helmholtz equation in a bounded domain with 
the appropriate boundary conditions. We have assumed here that periodic solutions exist. 
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A second assumption is that the solutions of the initial value problems with periodic bound- 
ary conditions approach to a limiting state as / — ► oo. While this is certainly expected for 
reasonably shaped scatterers, it should be proved in a precise treatment. 

The third assumption is the one we justified by appeal to computational experiments, which 
is that in the transverse magnetic (electric) case the electric (magnetic) held has no offset 
so that E(x, /) — e* lx ' t E(x) — ► 0 as t — ► oo (H(x, /) — e* lx ' t H(x) — ► 0) as t — ► oo. Since in the 
transverse magnetic (electric) problem E (H) is one-dimensional, this is really a question of 
the asymptotic behavior of the wave equation with periodic Dirichlet (Neumann) condition 
on the inner boundary and the absorbing boundary condition on the outer boundary. This 
has been known for a while. For three dimensions the question is more delicate, because 
the equations are coupled by divergence conditions, and also because both the electric and 
the magnetic held may have spurious components. 

In our paper [4], which is a companion to this one, we provide a rigorous statement of 
these assumptions and the appropriate results are proved correct for the two dimensional 
problems. 


7. CONCLUSIONS 

We have shown the origin of the spurious magnetic and electric offset helds which are 
encountered when time marching Maxwell equations for periodic solutions. Based on the 
equations governing these helds several remedies were discussed. We found that the best 
remedy is to modify the incident held, gradually allowing it to reach its desired form in the 
way specihed in section 4. This technique reliably reduces the magnetic offset to a negligible 
size and avoids the need for postprocessing of the results of the computation. 
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